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Abstract 

Hypersonic  vehicles  are  subjected  to  extreme  acoustic,  thermal  and  mechanical  loading 
with  strong  spatial  and  temporal  gradients  and  for  extended  periods  of  time.  Long  duration, 
3-D  simulations  of  non-linear  response  of  these  vehicles,  is  prohibitively  expensive  using 
available  Finite  Element  Methods  and  algorithms.  This  report  presents  recent  advances  of 
a  Generalized  Finite  Element  Method  (GFEM)  for  multiscale  non-linear  simulations.  This 
method  is  able  to  handle  complex  non-linear  problems  such  as  those  exhibiting  softening 
in  the  load-displacement  curve.  Cohesive  fracture  models  lead  to  this  class  of  non-linear 
behavior,  which  are  significantly  more  computationally  expensive  than  in  the  case  of  linear 
elastic  fracture  mechanics.  In  this  novel  GFEM,  scale -bridging  enrichment  functions  are 
updated  on  the  fly  during  the  non-linear  iterative  solution  process.  Non-linear  fine-scale 
solutions  are  embedded  in  the  global  scale  using  the  partition  of  unity  framework  of  the 
Generalized  FEM.  Damage  information  computed  at  fine-scale  problems  are  also  used  at 
the  coarse  scale  in  order  to  avoid  costly  non-linear  iterations  at  the  global  scale.  This 
method  enables  high-fidelity  non-linear  simulation  of  representative  aircraft  panels  using 
finite  element  meshes  that  are  orders  of  magnitude  coarser  than  those  required  by  available 
finite  element  methods. 

We  also  report  on  a  technique  to  perform  a  near-orthogonalization  of  scale-bridging  enrich¬ 
ments  used  in  the  multiscale  GFEM.  We  show  that,  for  any  discretization  error  level,  it  leads 
to  systems  of  equations  that  are  orders  of  magnitude  better  conditioned  than  in  available 
GFEMs.  This  so-called  Stable  Generalized  FEM  (SGFEM)  requires  minimal  modifications 
of  existing  GFEM  software  and  leads  to  optimal  convergence  rates,  regardless  of  the  pres¬ 
ence  of  singular  solutions  due  to  cracks.  We  also  show  that  the  error  within  enrichment 
zones  in  the  SGFEM  is  lower  than  in  the  GFEM.  This  is  important  for  fracture  mechanics 
problems  since  parameters  such  as  stress  intensity  factors  are  extracted  within  these  zones. 
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Generalized  Finite  Element  Approximations 


Generalized  FEM  approximation  spaces  (i.e.,  trial  and  test  spaces)  consist  of  three  compo¬ 
nents  -  (a)  patches  or  clouds,  (b)  a  partition  of  unity,  and  (c)  the  patch  or  cloud  approxima¬ 
tion  spaces.  The  main  ideas  of  the  GFEM  are  summarized  in  this  section. 

Consider  the  usual  finite  element  partitioning  flh  of  a  given  domain  Q,  in  which  flh  is  the 
union  of  individual  finite  elements  Qe.  e  —  1, . . .  ,nel.  The  basic  idea  behind  the  GFEM 
is  to  hierarchically  enrich  a  low-order  standard  FEM  approximation  space,  E>fem,  with 
special  functions  tailored  for  the  physics  of  the  problem  at  hand.  These  functions  belong 
to  a  space  E>enr  defined  using  the  partition  of  unity  property  of  Lagrangian  finite  element 
shape  functions,  i.e., 

X>a  =  i,  (i) 

aeJh 

where  a,  a  e  Jh  —  {1, . . .  ,uq},  is  the  index  of  a  node  in  a  finite  element  mesh  with  hq 
nodes.  Linear  FEM  shape  functions  are  adopted  in  this  work. 

The  test  and  trial  GFEM  space  f$>GFEM  are  given  by 

&GFEM  —  &FEM  +  §ENR  (2) 


where 


§fem  =  and  SENR=  £  NaXa ,  Xa  =  ttLaidai,  da,dai  e  M3  (3) 

aeJh  cceJ'e'n,.  ' 

Here,  i,  i  =  1 , . . . ,  n®nr,  denotes  the  index  of  the  enrichment  function  Lai  at  node  a  and  n®nr 
is  the  number  of  enrichments  at  node  a.  Enrichments  Lai .  i  =  1, . . .  .nfnr,  form  a  basis  of 
the  patch  or  cloud  space  Xa  ( )  with  ooa  being  the  support  of  the  FEM  shape  function  Na . 
The  set  Jhenr  C  Jh  has  the  indexes  of  the  nodes  with  enrichment  functions.  It  is  noted  that 
each  patch  (node)  may  adopt  a  different  basis,  depending  on  the  behavior  of  the  solution  of 
the  problem  over  the  node  support. 

Based  on  the  above  definitions,  a  GFEM  shape  function  at  a  node  a  G  Jhenr  is  given  by 

(f)ai(x)  =  Na(x)Lai(x)  (no  summation  on  a).  (4) 

The  definition  of  shape  functions  as  described  above  provides  great  flexibility  since  en¬ 
richment  functions  are  not  limited  to  polynomials  as  in  the  standard  FEM.  For  example, 
in  the  case  of  cohesive  fracture  problems  considered  in  this  study,  Heaviside  functions 
are  adopted  to  represent  discontinuities  arbitrarily  located  in  a  finite  element  mesh.  Fur¬ 
thermore,  enrichment  functions  for  multiscale  and  non-linear  problems  can  be  computed 
numerically  as  described  next. 
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Bridging  Scales  with  the  Generalized  Finite  Element  Method 


The  Generalized  Finite  Element  Method  with  global-local  enrichments  ( GFEMgl )  com¬ 
bines  ideas  from  the  classical  global-local  finite  element  method  with  the  Generalized  FEM 
described  in  the  previous  section.  In  contrast  to  available  Generalized  or  extended  FEMs, 
which  use  analytical  enrichment  functions,  this  method  provides  a  framework  that  allows 
the  enrichment  of  the  GFEM  solution  space  with  functions  obtained  from  the  solution  of 
local  boundary  value  problems.  The  boundary  conditions  for  local  problems  are  obtained 
from  the  solution  of  the  global  problem  discretized  with  a  coarse  finite  element  mesh. 
The  local  problems  can  be  accurately  solved  using  an  adaptive  GFEM ,  and  therefore  the 
GFEMgl  can  be  applied  to  problems  with  limited  a  priori  knowledge  about  the  solution 
like  those  involving  3-D  complex  fractures,  multiscale  or  non-linear  phenomena.  In  this 
method,  the  patch  or  cloud  approximation  spaces  are  built  with  the  aid  of  local  boundary 
value  problems  defined  in  a  neighborhood  EE  of  a  crack  or  other  local  feature  of  interest. 
Global-local  enrichment  functions  can  be  built  for  many  classes  of  problems.  Here  we 
report  on  a  three-dimensional  formulation  developed  for  propagating  non-linear  cohesive 
fractures.  Further  details  can  be  found  in  [12,  14].  It  is  noted  that  the  GFEM  developed 
in  project  FA9550-09-1-0401  assumed  linear  behavior  for  propagating  fractures  or  non¬ 
linear  but  stationary  fractures.  In  contrast,  the  GFEM  described  here  can  handle  the  case  of 
propagating  non-linear  fractures  with  load-displacement  curves  exhibiting  softening.  This 
creates  several  challenges  for  a  multiscale  method  since  it  requires  algorithms  able  to  deal 
with  load-dependent  discretization  spaces  while  avoiding  mapping  of  solutions  between 
spaces. 


Model  boundary  value  problem 


Fet  a  domain  Q.q  with  discontinuity  surfaces  Fcoh  be  occupied  by  a  body  to  be  open,  and 
bounded  by  a  smooth  boundary  Tq  that  involves  F^  and  f  [i  for  prescribed  displacement  ii 
and  traction  t,  respectively.  Figure 

The  body  can  be  characterized  by  a  single  variable,  the  displacement  field  uq  :  Qq  — *  Wldim 
(with  ridim  —  3  for  three  dimensions)  which  weakly  satisfies  equilibrium  in  a  Hilbert  space 
H1  as 


I  Or 


Vs  (8uq)  :  (7  {uq)  dV  +  5[»G]-fcoft(tt»G  1)  dS  +  T?  /  5»G-«GdS 


jjncoh 


ru 

7iG 


(5) 


=  /  duo  •  b  dV  +  /  8uc  tdS  +  rj  /  8ug  ■  u  dS 
JnG 


,lg 


ru 

7iG 


for  all  8uq  e  Hl.  Here,  we  use  notations  a  for  the  Cauchy  stress  tensor,  b  for  the 
volumetric  body  force,  p  for  the  penalty  parameter,  and  [i/g]  for  the  jump  of  displacement 
on  rcoh,  respectively. 

The  constitutive  relation  between  the  cohesive  traction,  tco/7,  and  the  displacement  jump, 
[hgJ,  is  in  general  highly  non-linear.  Global-local  enrichments  able  to  approximate  this 
behavior  are  presented  next. 
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Global  boundary  value  problem 


Local  boundary  value  problem 


Enrichments 


Figure  1:  GFEMgi  framework  for  two-scale  simulations  of  propagating  cohesive  fractures.  Numerically 
computed  solutions  of  the  extracted  local  boundary  value  problem  are  used  to  enrich  global  shape  functions 
while  solutions  of  the  original  global  boundary  value  problem  provide  boundary  conditions  for  the  local 
problem,  thus  defining  two  interdependent  problems  at  different  scales. 

Fine-scale  non-linear  local  problem 

Let  ukQl  G  E>q1(Qq)  be  a  GFEM  approximation  of  the  solution  uq  of  Problem  (5)  at 
the  (k  —  l)th  load/displacement  step.  Hereafter  a  load  and/or  displacement  step  is  denoted 
simply  by  load  step.  The  definitions  of  a  global  problem  to  compute  ukG  1  and  the  solution 
space  Sk-;  1  are  provided  later.  Global-local  enrichments  are  used  in  the  definition  of  1 . 
These  functions  are  the  solution  of  non-linear  local  problems  as  described  next. 

Let  a  sub-domain  Q.\  C  k2o  containing,  for  simplicity,  the  entire  pre-defined  crack  path  as 
illustrated  in  Figure  1.  Prescribed  displacements  uk  and  tractions  tk  at  the  kth  solution  step 
are  prescribed  on  TL  and  if  fli^,  respectively,  where  TL  denotes  the  boundary  of  Q\  . 
The  boundary  conditions  prescribed  on  Jl\  (TL  D  (T^  UT^))  are  provided  by  an  estimate, 
Uq  0,  of  the  global  solution  at  the  kth  load  step  and  defined  as 


k 
k  — 


,k- 

G  • 


(6) 


Solution  Uq  0  is  used  as  boundary  conditions  on  the  portion  of  Jl  that  does  not  intersect 
with  the  boundary  of  the  global  domain  Qc,-  This  is  a  key  aspect  of  the  method. 

Using  the  above  definitions,  the  weak  statement  of  the  non-linear  local  problem  at  the  kth 
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load  step  reads:  Find  ukL  e  S*  (T2l)  C  Hl  (T2l)  such  that  for  all  8ukL  e  (.Ql), 


Vs(5m|)  :  <t(b£)  dV  + 
+k 

+?7 


J  j~"  coh 


5[»il ■*«*(["£])  dS  +  T7 


8uk,  ■  uk,  dS 


/rLnrGu 


/rL\(rLn(rGuurcl)) 

4 


dS  = 


'A. 


dV  + 


Jin  12 


lL  '  UL 

8ukL-f  dS 


<5ttr • uk  dS  + 


/rLnrGu 


/jL\(rLn(i3fui3)) 


5h|  ■  \t(ukG  Q)  +  Kukc  o\  dS 


(7) 


for  the  penalty  parameter  77  and  the  spring  constant  K  defined  on  I  j  Di^1  and 
f L\  (Jl  H  (F^  UJq)),  respectively.  The  local  solution  space  (f^ij  is  defined  using  the 
standard  GFEM  shape  functions.  Much  finer  meshes  are  typically  used  at  local  than  in  the 
global  problem  as  illustrated  in  Figure  2.  This  figures  shows  the  application  of  the  GFEMgl 
to  a  three -point  bending  test.  This  problem  was  used  in  [12]  to  validate  the  proposed 
multiscale  method. 


Coarse-scale  global  problem 


Fine-scale  local  problem 


• :  nodes  enriched  by  polynomials  & 
Heaviside  functions 
• :  nodes  enriched  by  only  polynomials 


• :  nodes  enriched  by  only  polynomials 


Interface  between  local  domain 
boundary  and  global  domain 


parameters 


nodes  enriched  by  polynomials  & 
local  solutions 


Figure  2:  Model  problem  used  to  illustrate  the  non-linear  GFEMgi.  The  solution  computed  on  the  coarse 
global  mesh  provides  boundary  conditions  for  the  extracted  local  domain  in  a  neighborhood  of  non-linear 
cohesive  fracture.  A  fine  mesh  is  required  to  resolve  fine-scale  features  in  the  local  problem,  whereas  a 
coarse  mesh  is  used  to  capture  smooth  structural  behavior  in  the  global  problem.  Red  spheres  denote  nodes 
enriched  by  local  solutions  in  the  global  mesh  and  nodes  enriched  by  Heaviside  functions  in  the  local  mesh, 
respectively.  It  is  noted  that  the  local  mesh  does  not  match  the  global  one  at  the  local  domain  boundary. 


Global-Local  Enrichment  Functions  for  Cohesive  Fractures 

The  local  solution  ukL  defined  in  the  previous  section  is  used  to  build  the  following  GFEM 
shape  function  for  the  approximation  of  global  solution  uq  of  Problem  (5) 

($)a'k  =  NaukL,  (8) 
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where  the  partition  of  unity  function,  Na,  is  provided  by  a  coarse  global  FE  mesh  for 
Q(;  and  urL  has  the  role  of  an  enrichment  or  basis  function  for  the  patch  space  Xa(<Da)- 
Hereafter,  ukL  is  denoted  a  global- local  enrichment  function  at  the  kth  load  step.  The  cor¬ 
responding  global  GFEM  space  is  given  by  hierarchically  augmenting  the  standard  FEM 
solution  space  with  an  enrichment  space  ^nr  containing  shape  function  (j)a,k ,  i.e., 

—  ^fem  +  ^enr  =  ^fem  +  {NaUa,k  (no  summation  on  a),  a  e  J^sI}  (9) 

where  has  the  indexes  of  nodes  (patches)  enriched  with  global-local  functions.  A  node 
a  can  belong  to  set  J^gI  only  if  (Oa  C  fit.  Vector  u a'k  belongs  to  Xa(o)a)  and  is  given  by 

k  k,<u>  ) 

uau, 

44’<v>  \  do) 

k  k,<w> 

wauL  ) 

where  uf<u> ,  uR<v> ,  itf  u  are  the  components  of  the  local  solution  ukL  vector  in  the 
global  Cartesian  coordinate  directions  and  uka.  vkn .  wka  G  M  are  global  degrees  of  freedom. 
Equation  (10)  implies  that  G-L  enrichments  lead  to  only  three  additional  DOFs  per  global 
node,  regardless  of  the  size  of  the  local  problem  used  in  the  computation  of  urL. 

The  global  GFEM  space  fJf  defined  above  can  be  used  to  discretize  the  non-linear  global 
problem  (5)  and  find  a  global  approximation  ukG  e  SG ( Qq )  at  the  kth  load  step.  The  method¬ 
ology  is  illustrated  in  Figure  2.  The  global  solution  provides  boundary  conditions  for  fine- 
scale  problems  while  their  solutions  are  used  as  enrichment  functions  for  the  coarse-scale 
problem  through  the  partition  of  unity  framework  of  the  GFEM.  Figure  3  shows  the  load- 
displacement  curves  computed  with  the  GFEMgl  and  reference  numerical  and  experimental 
data  [12].  It  is  noted  that  the  GFEMgl  model  has  about  10  times  fewer  degrees  of  freedom 
than  in  the  case  of  available  adaptive  methods. 

It  is  noted  that  the  solution  space  SkG  is  adaptive:  It  changes  at  every  load  step  in  order 
to  approximate  the  non-linear  response  of  the  problem  while  keeping  the  global  mesh  un¬ 
changed.  This  change  must  be  properly  handled  when  solving  the  non-linear  equations  us¬ 
ing,  for  example,  Newton-Rhapson  algorithms.  In  particular,  the  vector  with  global  DOFs 
dk.  1  computed  at  the  previous  load  step  is  not  a  robust  choice  for  the  initialization  of  the 
Newton-Rhapson  non-linear  iterations  at  load  step  k:  The  global  vectors  dL(.  1  and  dG  repre¬ 
sent  coefficients  of  different  sets  of  GFEM  shape  functions.  Even  though  the  global  GFEM 
mesh  does  not  change,  the  solution  space  does.  An  efficient  and  robust  algorithm  to  deal 
with  this  issue  is  presented  in  [12].  It  is  based  on  the  solution  of  a  linear  problem  using 
a  secant  material  stiffness  instead  of  the  tangent  stiffness.  Further  details  can  be  found  in 
[12,  14]. 
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Figure  3:  Representative  problem  solved  with  the  non-linear  GFEMgl  [12]:  (normalized)  reaction  P  versus 
crack  mouth  opening  displacement  (CMOD)  curves  for  the  problem  shown  in  Figure  2.  The  GFEM gl  solution 
computed  in  the  global  problem  is  compared  with  available  reference  data.  FEM,  finite  element  method;  hp- 
GFEM,  hp  version  of  the  generalized  finite  element  method;  DOFs,  degrees  of  freedom;  GFEM gl,  generalized 
finite  element  method  with  globallocal  enrichments. 

Stable  Generalized  Finite  Element  Method 

In  this  section,  we  report  on  a  technique  to  perform  a  near-orthogonalization  of  enrichments 
used  in  the  Generalized  FEM.  The  technique  involves  modifications  to  enrichments  for  the 
GFEM  in  order  to  create  functions  that  are  near  orthogonal  to  the  finite  element  partition 
of  unity.  This  so-called  Stable  GFEM  ( SGFEM )  was  originally  proposed  by  Babuska  and 
Banerjee  [BB12].  They  have  shown  that  the  conditioning  of  the  SGFEM  is  not  worse  than 
that  of  the  standard  FEM.  In  this  project,  extensions  of  the  SGFEM  to  three-dimensional 
fracture  problems  were  developed  in  collaboration  with  Prof.  Ivo  Babuska  from  University 
of  Texas  at  Austin  and  Prof.  Uday  Banerjee  from  Syracuse  University.  This  collaboration 
is  at  no  cost  to  the  AFOSR.  A  summary  of  the  method  is  presented  below.  Details  on  this 
3-D  SGFEM  are  described  in  [10,  8]. 

In  the  SGFEM,  the  enrichment  functions  employed  in  GFEM  are  locally  modified  to  con¬ 
struct  the  patch  approximation  spaces  xa,  a  E  Jhenr.  The  modified  SGFEM  enrichment 
functions  La,(x)  E  Xa(o)a)  are  given  by 

Lai( x)  =  La\x)  -  lcoa  (Lai)  (x)  and  =  span{ Lai}<[  (11) 

where  l0ja  (Fai)  is  the  piecewise  tri-linear  finite  element  interpolant  of  the  enrichment  func¬ 
tion  Lai  on  the  patch  coa-  The  interpolant  I ma(La,)(%)  at  master  coordinate  £  of  a  finite 
element  T  with  nodes  J^(  t)  and  belonging  to  patch  (oa  is  given  by 

M£“')(4)=  E  d2) 
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where  vector  xp  has  the  coordinates  of  node  /3  of  element  t  and  N&  is  the  piecewise  tri- 
linear  FE  shape  function  for  node  j8.  Further  details  can  be  found  in  [10].  The  global 
enrichment  space  associated  with  Xa  is  denoted  by  §enr-  Therefore,  the  SGFEM  trial 
space  Ssgfem  is  given  by 

§  SGFEM  —  §  FEM  +  §ENR  ■  (13) 

The  SGFEM  shape  functions  (j)ai( x)  belonging  to  § enr  are  constructed  using  the  same 
framework  as  in  the  GFEM  and  are  given  by 


<p(jc)  =Na{x)Lai{x).  (14) 

Figure  4  illustrates  the  computation  of  SGFEM  enrichment  functions  and  shape  functions 
in  $enr  in  a  2-D  setting. 
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Figure  4:  Figure  illustrating  the  computation  of  an  SGFEM  enrichment  function  in  2-D.  The  picture  on  the 
left  shows  the  construction  of  a  GFEM  shape  function.  The  center  picture  features  the  original  enrichment 
function,  Lai,  at  the  top,  the  piecewise  bi-linear  finite  element  interpolant  of  which  is  in  the  middle,  Iffla  (Lai), 
and  the  modified  SGFEM  enrichment  function,  La\  is  shown  at  the  bottom.  The  picture  on  the  right  shows 
the  construction  of  an  SGFEM  shape  function,  0  . 
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